scalar = "dti_md"
HCPD_ageeffects <- read.csv(paste0(output_root, "/", "HCPD", "/GAM/", scalar, "/HCPD_GAM_dev_measures.csv"))
HBN_ageeffects <- read.csv(paste0(output_root, "/", "HBN", "/GAM/", scalar, "/HBN_GAM_dev_measures.csv"))
PNC_ageeffects <- read.csv(paste0(output_root, "/", "PNC", "/GAM/", scalar, "/PNC_GAM_dev_measures.csv"))# make a list of the different df names: [dataset]_GAM_ageeffects_[scalar]
datasets <- c("HCPD", "HBN", "PNC")
ageeffect_df_names <- c()
for (dataset in datasets) {
ageeffect_df_names <- append(ageeffect_df_names, paste0(dataset, "_ageeffects"))
}
# format age effect df's
ageeffect_dfs <- lapply(ageeffect_df_names, format_ageeffect)## There are 2392/2400 significant nodes
## There are 2350/2400 significant nodes
## There are 2273/2400 significant nodes
names(ageeffect.fdr_dfs) <- ageeffect_df_names
tracts <- setdiff(unique(ageeffect.fdr_dfs$HCPD_ageeffects$tract_label), "Corticospinal") # CST for supplementcor_HCPD_HBN <- round(cor.test(ageeffect.fdr_dfs$HCPD_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$HBN_ageeffects$GAM.smooth.AdjRsq)$estimate, 2)
cor_HCPD_PNC <- round(cor.test(ageeffect.fdr_dfs$HCPD_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$PNC_ageeffects$GAM.smooth.AdjRsq)$estimate, 2)
cor_HBN_PNC <- round(cor.test(ageeffect.fdr_dfs$HBN_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$PNC_ageeffects$GAM.smooth.AdjRsq)$estimate, 2)
# permute for p-values
HCPD_HBN_ageeffect_perm <- perm_cor_test(ageeffect.fdr_dfs$HCPD_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$HBN_ageeffects$GAM.smooth.AdjRsq)
HCPD_PNC_ageeffect_perm <- perm_cor_test(ageeffect.fdr_dfs$HCPD_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$PNC_ageeffects$GAM.smooth.AdjRsq)
PNC_HBN_ageeffect_perm <- perm_cor_test(ageeffect.fdr_dfs$PNC_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$HBN_ageeffects$GAM.smooth.AdjRsq)
HCPD_to_corr <- ageeffect.fdr_dfs$HCPD_ageeffects %>% mutate(HCPD = abs(GAM.smooth.AdjRsq)) %>% mutate(tract_node = paste0(tractID, "_", nodeID)) %>% select(tract_node, HCPD)
HBN_to_corr <- ageeffect.fdr_dfs$HBN_ageeffects %>% mutate(HBN = abs(GAM.smooth.AdjRsq)) %>% mutate(tract_node = paste0(tractID, "_", nodeID)) %>% select(tract_node, HBN)
PNC_to_corr <- ageeffect.fdr_dfs$PNC_ageeffects %>% mutate(PNC = abs(GAM.smooth.AdjRsq)) %>% mutate(tract_node = paste0(tractID, "_", nodeID)) %>% select(tract_node, PNC)
HCPD_HBN_ageeffects <- merge(HCPD_to_corr, HBN_to_corr, by = "tract_node")
HCPD_PNC_ageeffects <- merge(HCPD_to_corr, PNC_to_corr, by = "tract_node")
HBN_PNC_ageeffects <- merge(HBN_to_corr, PNC_to_corr, by = "tract_node")
text_HCPD_HBN <- bquote(italic(r) == .(cor_HCPD_HBN) * "," ~ italic(p[perm]) < .(format(round(HCPD_HBN_ageeffect_perm$p_value, 4), scientific = FALSE)))
HCPD_HBN_ageeffects_plot <- hex_plot(HCPD_HBN_ageeffects, "HCPD", "HBN", text_HCPD_HBN, ylim1 = 0 , ylim2 = 0.64, xlim1 = 0, xlim2 = 0.62, x_text = 0, y_text = 0.64)
text_HCPD_PNC <- bquote(italic(r) == .(cor_HCPD_PNC) * "," ~ italic(p[perm]) < .(format(round(HCPD_PNC_ageeffect_perm$p_value, 4), scientific = FALSE)))
HCPD_PNC_ageeffects_plot <- hex_plot(HCPD_PNC_ageeffects, "HCPD", "PNC", text_HCPD_PNC, ylim1 = 0 , ylim2 = 0.64, xlim1 = 0, xlim2 = 0.62, x_text = 0, y_text = 0.64)
text_HBN_PNC <- bquote(italic(r) == .(cor_HBN_PNC) * "," ~ italic(p[perm]) < .(format(round(PNC_HBN_ageeffect_perm$p_value, 4), scientific = FALSE)))
HBN_PNC_ageeffects_plot <- hex_plot(HBN_PNC_ageeffects, "HBN", "PNC", text_HBN_PNC, ylim1 = 0 , ylim2 = 0.64, xlim1 = 0, xlim2 = 0.62, x_text = 0, y_text = 0.64)
# make bin colors the same across plots
# get the maximum bin counts from all datasets
max_count_HCPD_HBN <- max(ggplot_build(HCPD_HBN_ageeffects_plot)$data[[1]]$count, na.rm = TRUE)
max_count_HCPD_PNC <- max(ggplot_build(HCPD_PNC_ageeffects_plot)$data[[1]]$count, na.rm = TRUE)
max_count_HBN_PNC <- max(ggplot_build(HBN_PNC_ageeffects_plot)$data[[1]]$count, na.rm = TRUE)
global_max_count <- max(max_count_HCPD_HBN, max_count_HCPD_PNC, max_count_HBN_PNC)
scale_limits <- c(0, global_max_count)
HCPD_HBN_ageeffects_plot <- HCPD_HBN_ageeffects_plot +
paletteer::scale_fill_paletteer_c("ggthemes::Blue-Green Sequential",
direction = -1,
limits = scale_limits,
oob = scales::squish)
HCPD_PNC_ageeffects_plot <- HCPD_PNC_ageeffects_plot +
paletteer::scale_fill_paletteer_c("ggthemes::Blue-Green Sequential",
direction = -1,
limits = scale_limits,
oob = scales::squish)
HBN_PNC_ageeffects_plot <- HBN_PNC_ageeffects_plot +
paletteer::scale_fill_paletteer_c("ggthemes::Blue-Green Sequential",
direction = -1,
limits = scale_limits,
oob = scales::squish)
ageeffects_plots <- ggarrange(HCPD_PNC_ageeffects_plot, HBN_PNC_ageeffects_plot, HCPD_HBN_ageeffects_plot, nrow = 1, common.legend=T, legend = c("bottom"), labels = "auto")
title.grob <- textGrob("Age Effect Magnitude of Mean Diffusivity Across Datasets",
gp=gpar(col="black", fontsize = 24, face = "bold"), hjust = 0.5)
grid.arrange(arrangeGrob(ageeffects_plots, top = title.grob))PNC <- readRDS(sprintf("%1$s/%2$s/tract_profiles/tract_profiles_for_NEST.RData", output_root, "PNC"))
HCPD <- readRDS(sprintf("%1$s/%2$s/tract_profiles/tract_profiles_for_NEST.RData", output_root, "HCPD"))
HBN <- readRDS(sprintf("%1$s/%2$s/tract_profiles/tract_profiles_for_NEST.RData", output_root, "HBN"))tract profiles coefficient of variance
datasets <- list(HCPD = HCPD, HBN = HBN, PNC = PNC)
# create summary dfs
for (dataset_name in names(datasets)) {
dataset <- datasets[[dataset_name]]
make_summary_dfs("dti_md", dataset_name, dataset) # this creates summary_[dataset]_[scalar] dataframes that include mean, sd, se
}clipEnds = 5
bin_size = 5
PNC_cv <- summary_PNC_dti_md %>%
filter((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) |
(nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size)) |
(nodeID > (50-bin_size) & nodeID <= (50+bin_size))) %>%
mutate(node_position = case_when((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) |
nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size) ~ "superficial",
(nodeID > (50-bin_size) & nodeID <= (50+bin_size)) ~ "deep")) %>%
group_by(tract_label, tractID, node_position) %>%
summarise(average_cv = round(mean(cv),2)) %>%
pivot_wider(names_from = node_position, values_from = average_cv, names_prefix = "average_cv_") %>%
select(tractID, average_cv_deep, average_cv_superficial) %>%
rename(`PNC Deep` = average_cv_deep,
`PNC Superficial` = average_cv_superficial,
Tract = tractID)
PNC_cv$Tract <- gsub("_", " ", PNC_cv$Tract)
PNC_cv <- PNC_cv[,-1]
PNC_cv <- PNC_cv %>% ungroup()
HCPD_cv <- summary_HCPD_dti_md %>%
filter((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) |
(nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size)) |
(nodeID > (50-bin_size) & nodeID <= (50+bin_size))) %>%
mutate(node_position = case_when((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) |
nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size) ~ "superficial",
(nodeID > (50-bin_size) & nodeID <= (50+bin_size)) ~ "deep")) %>%
group_by(tract_label, tractID, node_position) %>%
summarise(average_cv = round(mean(cv),2)) %>%
pivot_wider(names_from = node_position, values_from = average_cv, names_prefix = "average_cv_") %>%
select(tractID, average_cv_deep, average_cv_superficial) %>%
rename(`HCPD Deep` = average_cv_deep,
`HCPD Superficial` = average_cv_superficial,
Tract = tractID)
HCPD_cv$Tract <- gsub("_", " ", HCPD_cv$Tract)
HCPD_cv <- HCPD_cv[,-1]
HCPD_cv <- HCPD_cv %>% ungroup()
HBN_cv <- summary_HBN_dti_md %>%
filter((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) |
(nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size)) |
(nodeID > (50-bin_size) & nodeID <= (50+bin_size))) %>%
mutate(node_position = case_when((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) |
nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size) ~ "superficial",
(nodeID > (50-bin_size) & nodeID <= (50+bin_size)) ~ "deep")) %>%
group_by(tract_label, tractID, node_position) %>%
summarise(average_cv = round(mean(cv),2)) %>%
pivot_wider(names_from = node_position, values_from = average_cv, names_prefix = "average_cv_") %>%
select(tractID, average_cv_deep, average_cv_superficial) %>%
rename(`HBN Deep` = average_cv_deep,
`HBN Superficial` = average_cv_superficial,
Tract = tractID)
HBN_cv$Tract <- gsub("_", " ", HBN_cv$Tract)
HBN_cv <- HBN_cv[,-1]
HBN_cv <- HBN_cv %>% ungroup()
combined_cv <- left_join(PNC_cv, HCPD_cv)
combined_cv <- left_join(combined_cv, HBN_cv)
colnames(combined_cv) <- sub("^[^ ]+ ", "", colnames(combined_cv))
combined_cv$Tract <- gsub("Left", "L", combined_cv$Tract)
combined_cv$Tract <- gsub("Right", "R", combined_cv$Tract)
combined_cv <- rbind(combined_cv, c("Average Across Tracts", round(colMeans(combined_cv[,c(2:7)]), 2)))
cv_table_html <- combined_cv %>%
kbl(format = "html", align = "lcccccc", digits = 2) %>%
add_header_above(c(" " = 1, "PNC" = 2, "HCP-D" = 2, "HBN" = 2)) %>%
add_header_above(c(" " = 1, "Average Coefficient of Variation" = 6), bold = TRUE) %>%
kable_styling(full_width = FALSE, position = "center") %>%
kable_classic(full_width = FALSE, html_font = "Arial") %>%
row_spec(0:nrow(combined_cv), extra_css = "line-height: 0.5;") %>%
column_spec(c(2, 4, 6), background = "gray90", include_thead = TRUE)
cv_table_html| Tract | Deep | Superficial | Deep | Superficial | Deep | Superficial |
|---|---|---|---|---|---|---|
| L Arcuate | 4.26 | 4.11 | 3.8 | 3.87 | 7.83 | 5.25 |
| R Arcuate | 4.99 | 4.15 | 4.15 | 3.79 | 7.26 | 5.48 |
| Callosum Anterior Frontal | 5.06 | 3.53 | 4.41 | 3.54 | 8.29 | 4.83 |
| Callosum Motor | 11.42 | 4.63 | 8.34 | 3.77 | 10.97 | 7.04 |
| Callosum Occipital | 7.36 | 4.34 | 4.43 | 4.25 | 14.79 | 5.7 |
| Callosum Orbital | 6.98 | 4.08 | 4.23 | 3.87 | 12.04 | 6.37 |
| Callosum Posterior Parietal | 5.56 | 4.87 | 3.9 | 3.74 | 12.97 | 6.28 |
| Callosum Superior Frontal | 10.38 | 4.76 | 6.69 | 3.84 | 8.8 | 6.25 |
| Callosum Superior Parietal | 8.13 | 4.23 | 4.52 | 3.48 | 11.65 | 6.32 |
| Callosum Temporal | 8.61 | 4.76 | 5.53 | 3.82 | 14.68 | 9.45 |
| L Corticospinal | 3.1 | 10.84 | 3.22 | 6.24 | 7.33 | 8.74 |
| R Corticospinal | 3.4 | 10.73 | 2.87 | 6.56 | 7.02 | 8.86 |
| L Inferior Fronto.occipital | 9.67 | 4.58 | 6.01 | 4.2 | 7.59 | 5.48 |
| R Inferior Fronto.occipital | 5.55 | 4.95 | 5.32 | 4.37 | 7.16 | 5.49 |
| L Inferior Longitudinal | 4.79 | 4.66 | 4.09 | 3.93 | 6.5 | 6.43 |
| R Inferior Longitudinal | 5.7 | 4.57 | 4.14 | 3.81 | 6.87 | 5.62 |
| L Posterior Arcuate | 4.97 | 4.21 | 4.22 | 3.66 | 6.1 | 4.37 |
| R Posterior Arcuate | 5.19 | 4.59 | 4.1 | 3.65 | 6.22 | 4.9 |
| L Superior Longitudinal | 4.25 | 4.35 | 3.78 | 3.89 | 7.56 | 5.22 |
| R Superior Longitudinal | 4.49 | 4.35 | 3.46 | 3.81 | 7.35 | 5.43 |
| L Uncinate | 3.45 | 4.51 | 4.38 | 4.13 | 7.2 | 6.91 |
| R Uncinate | 3.29 | 4.27 | 3.44 | 4.14 | 6.62 | 6.37 |
| L Vertical Occipital | 4.77 | 3.9 | 3.88 | 4.26 | 5.48 | 5.05 |
| R Vertical Occipital | 5.65 | 4.03 | 3.94 | 4.05 | 5.57 | 5.18 |
| Average Across Tracts | 5.88 | 4.92 | 4.45 | 4.11 | 8.49 | 6.13 |
cst_PNC <- readRDS("/cbica/projects/luo_wm_dev/two_axes_manuscript/output/PNC/NEST/dti_md/Corticospinal_bin5_clip5_1sided.RData")
NEST_PNC <- data.frame(c("Corticospinal", cst_PNC$pval, "PNC")) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))
cst_HCPD <- readRDS("/cbica/projects/luo_wm_dev/two_axes_manuscript/output/HCPD/NEST/dti_md/Corticospinal_bin5_clip5_1sided.RData")
NEST_HCPD <- data.frame(c("Corticospinal", cst_HCPD$pval, "HCP-D")) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))
cst_HBN <- readRDS("/cbica/projects/luo_wm_dev/two_axes_manuscript/output/HBN/NEST/dti_md/Corticospinal_bin5_clip5_1sided.RData")
NEST_HBN <- data.frame(c("Corticospinal", cst_HBN$pval, "HBN")) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))ageeffect_measure = "GAM.smooth.AdjRsq"
colors <- c("#3FB8AFFF", "#FF9E9DFF", "#cc0468")
names(colors) <- c("HCP-D", "HBN", "PNC")
bin_num_nodes = 5
clipEnds = 5
# prepare the data for making lollipop plots
df_all <- bind_rows(
ageeffect.fdr_dfs$HCPD_ageeffects %>% mutate(Dataset = "HCP-D"),
ageeffect.fdr_dfs$HBN_ageeffects %>% mutate(Dataset = "HBN"),
ageeffect.fdr_dfs$PNC_ageeffects %>% mutate(Dataset = "PNC")
)
df_all$Dataset <- factor(df_all$Dataset, levels = c("HCP-D", "PNC", "HBN"))
df_formatted <- format_deep_superficial(df_all, ageeffect_measure, bin_num_nodes, clipEnds)
df_formatted <- df_formatted[complete.cases(df_formatted), ]
df_formatted <- df_formatted %>% filter(tract_label == "Corticospinal")
# merge with NEST results
NEST <- rbind(NEST_HCPD, NEST_HBN, NEST_PNC)
lollipop_data <- prepare_lollipop_data(df_formatted)
lollipop_data$Dataset <- factor(lollipop_data$Dataset, levels = c("PNC", "HCP-D", "HBN"))
# plot lollipops
lollipop_plot_CST <- make_lollipop_plot(tract = "Corticospinal", data = lollipop_data)ageeffect_measure = "GAM.smooth.AdjRsq"
ageeffects_plot_CST <- plot_age_effect(tract = "Corticospinal", scalar = "dti_md", ageeffect_measure = ageeffect_measure,
color_HCPD = color_HCPD, color_HBN = color_HBN, color_PNC = color_PNC,
clipEnds = 5, ylim2 = 0.41, ylim1 = 0, fontsize = 18, legend_position = "none")CST_plot <- plot_grid(ggdraw() + draw_label("Corticospinal", size = 18, y = 0.2, hjust = 0.5),
plot_grid(plotlist = list(ageeffects_plot_CST, lollipop_plot_CST),
align = "h", ncol = 2, rel_widths = c(1, 0.4)), ncol = 1, rel_heights = c(0.25, 1)) + bgcolor("white")
CST_plotscalar = "dti_fa"
PNC_ageeffects <- read.csv(paste0(output_root, "/", "PNC", "/GAM/", scalar, "/PNC_GAM_dev_measures.csv"))
HCPD_ageeffects <- read.csv(paste0(output_root, "/", "HCPD", "/GAM/", scalar, "/HCPD_GAM_dev_measures.csv"))
HBN_ageeffects <- read.csv(paste0(output_root, "/", "HBN", "/GAM/", scalar, "/HBN_GAM_dev_measures.csv"))# make a list of the different df names: [dataset]_GAM_ageeffects_[scalar]
datasets <- c("PNC", "HCPD", "HBN")
ageeffect_df_names <- c()
for (dataset in datasets) {
ageeffect_df_names <- append(ageeffect_df_names, paste0(dataset, "_ageeffects"))
}
# format age effect df's
ageeffect_dfs <- lapply(ageeffect_df_names, format_ageeffect)## There are 1621/2400 significant nodes
## There are 1061/2400 significant nodes
## There are 1947/2400 significant nodes
names(ageeffect.fdr_dfs) <- ageeffect_df_names
tracts <- setdiff(unique(ageeffect.fdr_dfs$HCPD_ageeffects$tract_label), "Corticospinal") # CST for supplement
# PNC: There are 1621/2400 significant nodes - 67.5%
# HCP-D: There are 1061/2400 significant nodes - 44.2%
# HBN: There are 1947/2400 significant nodes - 81.1%ageeffect_measure = "GAM.smooth.AdjRsq"
ageeffects_plots <- lapply(tracts, plot_age_effect, scalar = "dti_fa", ageeffect_measure = ageeffect_measure,
color_HCPD = color_HCPD, color_HBN = color_HBN, color_PNC = color_PNC, clipEnds = 5, ylim2 = 0.25, ylim1 = 0, fontsize
= 20, legend_position = "none")
names(ageeffects_plots) <- tractsarrange callosal plots
#ggsave(paste0(suppfig_dir, "supp_ageeffects_fig1_FA.svg"), arrange_callosum_plots_nolollipop(ageeffects_plots), height = 10, width = 20, units = "in")
#ggsave(paste0(suppfig_dir, "supp_ageeffects_fig1_FA.png"), arrange_callosum_plots_nolollipop(ageeffects_plots), height = 10, width = 20, units = "in")arrange association plots
ageeffect_measure = "GAM.smooth.AdjRsq"
ageeffects_plots <- lapply(tracts, plot_age_effect, scalar = "dti_fa", ageeffect_measure = ageeffect_measure,
color_HCPD = color_HCPD, color_HBN = color_HBN, color_PNC = color_PNC,
clipEnds = 5, ylim2 = 0.25, ylim1 = 0, fontsize = 18, legend_position = "none")
names(ageeffects_plots) <- tracts#ggsave(paste0(suppfig_dir, "supp_ageeffects_fig2_FA.svg"), arrange_association_plots_nolollipop(ageeffects_plots), height = 10, width = 20, units = "in")
#ggsave(paste0(suppfig_dir, "supp_ageeffects_fig2_FA.png"), arrange_association_plots_nolollipop(ageeffects_plots), height = 10, width = 20, units = "in")# reload dti_md files lol
scalar = "dti_md"
HCPD_ageeffects <- read.csv(paste0(output_root, "/", "HCPD", "/GAM/", scalar, "/HCPD_GAM_dev_measures.csv"))
HBN_ageeffects <- read.csv(paste0(output_root, "/", "HBN", "/GAM/", scalar, "/HBN_GAM_dev_measures.csv"))
PNC_ageeffects <- read.csv(paste0(output_root, "/", "PNC", "/GAM/", scalar, "/PNC_GAM_dev_measures.csv"))# make a list of the different df names: [dataset]_GAM_ageeffects_[scalar]
datasets <- c("HCPD", "HBN", "PNC")
ageeffect_df_names <- c()
for (dataset in datasets) {
ageeffect_df_names <- append(ageeffect_df_names, paste0(dataset, "_ageeffects"))
}
# format age effect df's
ageeffect_dfs <- lapply(ageeffect_df_names, format_ageeffect)## There are 2392/2400 significant nodes
## There are 2350/2400 significant nodes
## There are 2273/2400 significant nodes
names(ageeffect.fdr_dfs) <- ageeffect_df_names
tracts <- setdiff(unique(ageeffect.fdr_dfs$HCPD_ageeffects$tract_label), "Corticospinal") # CST for supplement# load glasser labels
glasser_labels <- read.csv("/cbica/projects/luo_wm_dev/atlases/glasser/HCP-MMP1_UniqueRegionList.csv")
glasser_labels$regionID <- c(1:360)
glasser_labels$region <- gsub("7Pl", "7PL", glasser_labels$region)
# load maps for depth = 1.5mm (disregard a,b, and c. Variables are just there to prevent lapply output from being printed lol)
a <- lapply("1.5", load_maps, "PNC") # makes lh_maps_PNC, rh_maps_PNC
b <- lapply("1.5", load_maps, "HCPD") # makes lh_maps_HCPD, rh_maps_HCPD
c <- lapply("1.5", load_maps, "HBN") # makes lh_maps_HBN, rh_maps_HBNglasser_SAaxis <- read.csv("/cbica/projects/luo_wm_dev/SAaxis/glasser_SAaxis.csv")
glasser_SAaxis <- glasser_SAaxis %>% select(SA.axis_rank, label)
glasser_SAaxis$regionName <- gsub("_ROI", "", glasser_SAaxis$label)
glasser_SAaxis$regionName <- gsub("^(.)_(.*)$", "\\2_\\1", glasser_SAaxis$region)threshold=0.3
PNC_deveffects_5_agemat <- ageeffect.fdr_dfs$PNC_ageeffects %>% filter((nodeID < 10 & nodeID > 4) | (nodeID < 95 & nodeID > 89)) %>% mutate(node_position = case_when((nodeID < 10 & nodeID > 4) ~ "end1", (nodeID < 95 & nodeID > 89) ~ "end2")) %>%
select(tract_label, tractID, nodeID, node_position, hemi, smooth.peak.change, smooth.decrease.offset, smooth.last.change, smooth.slowing.onset) %>%
group_by(tractID, node_position, hemi) %>% summarise(mean_peak_change = mean(smooth.peak.change, na.rm = T),
mean_ageeffect = mean(smooth.decrease.offset, na.rm = T),
mean_last_change = mean(smooth.last.change, na.rm = T),
mean_dev_slowing = mean(smooth.slowing.onset, na.rm = T))
HCPD_deveffects_5_agemat <- ageeffect.fdr_dfs$HCPD_ageeffects %>% filter((nodeID < 10 & nodeID > 4) | (nodeID < 95 & nodeID > 89)) %>% mutate(node_position = case_when((nodeID < 11 & nodeID > 4) ~ "end1", (nodeID < 95 & nodeID > 89) ~ "end2")) %>%
select(tract_label, tractID, nodeID, node_position, hemi, smooth.peak.change, smooth.decrease.offset, smooth.last.change, smooth.slowing.onset) %>%
group_by(tractID, node_position, hemi) %>% summarise(mean_peak_change = mean(smooth.peak.change, na.rm = T),
mean_ageeffect = mean(smooth.decrease.offset, na.rm = T), # mean_ageeffect = mean_age_mat
mean_last_change = mean(smooth.last.change, na.rm = T),
mean_dev_slowing = mean(smooth.slowing.onset, na.rm = T))
HBN_deveffects_5_agemat <- ageeffect.fdr_dfs$HBN_ageeffects %>% filter((nodeID < 10 & nodeID > 4) | (nodeID < 95 & nodeID > 89)) %>% mutate(node_position = case_when((nodeID < 10 & nodeID > 4) ~ "end1", (nodeID < 95 & nodeID > 89) ~ "end2")) %>%
select(tract_label, tractID, nodeID, node_position, hemi, smooth.peak.change, smooth.decrease.offset, smooth.last.change, smooth.slowing.onset) %>%
group_by(tractID, node_position, hemi) %>% summarise(mean_peak_change = mean(smooth.peak.change, na.rm = T),
mean_ageeffect = mean(smooth.decrease.offset, na.rm = T),
mean_last_change = mean(smooth.last.change, na.rm = T),
mean_dev_slowing = mean(smooth.slowing.onset, na.rm = T))
make_maps("PNC", "5_agemat") # makes [bundle_name]_deveffect_[dataset] for each dataset and bundle. e.g. IFOL_deveffect_PNC
make_maps("HCPD", "5_agemat")
make_maps("HBN", "5_agemat")region_all_PNC <- aggregate_age_maps("PNC")
lh_by_region_PNC <- region_all_PNC[[1]]
rh_by_region_PNC <- region_all_PNC[[2]]
region_all_HCPD <- aggregate_age_maps("HCPD")
lh_by_region_HCPD <- region_all_HCPD[[1]]
rh_by_region_HCPD <- region_all_HCPD[[2]]
region_all_HBN <- aggregate_age_maps("HBN")
lh_by_region_HBN <- region_all_HBN[[1]]
rh_by_region_HBN <- region_all_HBN[[2]]all_endpoints_PNC <- compute_mean_SA("PNC")
lh_all_endpoints_PNC <- all_endpoints_PNC[[1]]
rh_all_endpoints_PNC <- all_endpoints_PNC[[2]]
all_endpoints_PNC <- all_endpoints_PNC[[3]]
all_endpoints_HCPD <- compute_mean_SA("HCPD")
lh_all_endpoints_HCPD <- all_endpoints_HCPD[[1]]
rh_all_endpoints_HCPD <- all_endpoints_HCPD[[2]]
all_endpoints_HCPD <- all_endpoints_HCPD[[3]]
all_endpoints_HBN <- compute_mean_SA("HBN")
lh_all_endpoints_HBN <- all_endpoints_HBN[[1]]
rh_all_endpoints_HBN <- all_endpoints_HBN[[2]]
all_endpoints_HBN <- all_endpoints_HBN[[3]]
combined_df <- bind_rows(all_endpoints_PNC %>% mutate(dataset = "PNC"),
all_endpoints_HCPD %>% mutate(dataset = "HCPD"),
all_endpoints_HBN %>% mutate(dataset = "HBN"))
# Calculate the averages of mean_SA and age_effect across all dataframes
all_endpoints_avgdatasets <- combined_df %>%
group_by(bundle_name, end) %>% summarize(mean_SA = mean(mean_SA, na.rm = TRUE),
age_effect = mean(age_effect, na.rm = TRUE),
.groups = "drop")PNC_tractlevel_SA_pearsons <- read.csv('/cbica/projects/luo_wm_dev/two_axes_manuscript/output/PNC/suppfigs_spintests/supp_tractlevel_pearson_pspin_allendpoints.csv')
HCPD_tractlevel_SA_pearsons <- read.csv('/cbica/projects/luo_wm_dev/two_axes_manuscript/output/HCPD/suppfigs_spintests/supp_tractlevel_pearson_pspin_allendpoints.csv')
HBN_tractlevel_SA_pearsons <- read.csv('/cbica/projects/luo_wm_dev/two_axes_manuscript/output/HBN/suppfigs_spintests/supp_tractlevel_pearson_pspin_allendpoints.csv')
avgdatasets_tractlevel_SA_pearsons <- read.csv('/cbica/projects/luo_wm_dev/two_axes_manuscript/output/avgdatasets/suppfigs_spintests/supp_tractlevel_pearson_pspin_allendpoints.csv')
SA_age_mat_PNC <- plot_meanSA_by_age_mat("PNC", bquote(paste(italic("r"), " = ", .(round(PNC_tractlevel_SA_pearsons$rho.emp, 2)), ", ", italic("p"[spin]), " < 0.0001"))) # pval=0
SA_age_mat_HCPD <- plot_meanSA_by_age_mat("HCPD", bquote(paste(italic("r"), " = ", .(round(HCPD_tractlevel_SA_pearsons$rho.emp, 2)), ", ", italic("p"[spin]), " < 0.0001"))) # pval=0
SA_age_mat_HBN <- plot_meanSA_by_age_mat("HBN", bquote(paste(italic("r"), " = ", .(round(HBN_tractlevel_SA_pearsons$rho.emp, 2)), ", ", italic("p"[spin]), " = N.S.")))
SA_age_mat_avgdatasets <- plot_meanSA_by_age_mat("avgdatasets", bquote(paste(italic("r"), " = ", .(round(avgdatasets_tractlevel_SA_pearsons$rho.emp, 2)), ", ", italic("p"[spin]), " < 0.0001"))) + ggtitle("Average Across Datasets")
SA_age_mat_plot_final <- ggarrange(SA_age_mat_PNC, SA_age_mat_HCPD, SA_age_mat_HBN, SA_age_mat_avgdatasets, ncol = 4, common.legend = T, legend = "bottom")
x.grob <- textGrob("Mean S-A Axis Rank of Cortical Endpoint",
gp=gpar(col="black", fontsize=24))
y.grob <- textGrob("Age of Maturation (Years)",
gp=gpar(col="black", fontsize=24), rot=90)
grid.arrange(arrangeGrob(SA_age_mat_plot_final, left = y.grob, bottom = x.grob))#ggsave(paste0(suppfig_dir, "supp_SA_tract_level_pearson_allendpoints.png"), grid.arrange(arrangeGrob(SA_age_mat_plot_final, left = y.grob, bottom = x.grob)), height = 6, width = 22, units = "in")all_endpoints_PNC_binary <- all_endpoints_PNC
max_value <- max(all_endpoints_PNC_binary$age_effect, na.rm = TRUE)
all_endpoints_PNC_binary$age_effect[all_endpoints_PNC_binary$age_effect == max_value] <- NA
all_endpoints_PNC_binary <- all_endpoints_PNC_binary %>% mutate(maturation_status = ifelse(is.na(age_effect), 0, 1))
all_endpoints_HCPD_binary <- all_endpoints_HCPD
max_value <- max(all_endpoints_HCPD_binary$age_effect, na.rm = TRUE)
all_endpoints_HCPD_binary$age_effect[all_endpoints_HCPD_binary$age_effect == max_value] <- NA
all_endpoints_HCPD_binary <- all_endpoints_HCPD_binary %>% mutate(maturation_status = ifelse(is.na(age_effect), 0, 1)) # 0 = not matured, 1 = matured
all_endpoints_HBN_binary <- all_endpoints_HBN
max_value <- max(all_endpoints_HBN_binary$age_effect, na.rm = TRUE)
all_endpoints_HBN_binary$age_effect[all_endpoints_HBN_binary$age_effect == max_value] <- NA
all_endpoints_HBN_binary <- all_endpoints_HBN_binary %>% mutate(maturation_status = ifelse(is.na(age_effect), 0, 1))
all_endpoints_avgdatasets_binary <- all_endpoints_avgdatasets
max_value <- max(all_endpoints_avgdatasets_binary$age_effect, na.rm = TRUE)
all_endpoints_avgdatasets_binary$age_effect[all_endpoints_avgdatasets_binary$age_effect == max_value] <- NA
all_endpoints_avgdatasets_binary <- all_endpoints_avgdatasets_binary %>% mutate(maturation_status = ifelse(is.na(age_effect), 0, 1))
PNC_tractlevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/suppfigs_spintests/supp_tractlevel_pearson_ttest_pspin_maturedonly.csv")
HCPD_tractlevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/suppfigs_spintests/supp_tractlevel_pearson_ttest_pspin_maturedonly.csv")
HBN_tractlevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/suppfigs_spintests/supp_tractlevel_pearson_ttest_pspin_maturedonly.csv")
avgdatasets_tractlevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/avgdatasets/tract_profiles/suppfigs_spintests/supp_tractlevel_pearson_ttest_pspin_maturedonly.csv")
SA_age_mat_binary_PNC <- plot_meanSA_by_age_mat("PNC", bquote(paste(italic("r"), " = ", .(round(PNC_tractlevel_SA_ttest_p$rho.emp, 2)), ", ",
italic("p"[spin]), " < 0.0001")), binary = TRUE) # pval=0
SA_age_mat_binary_HCPD <- plot_meanSA_by_age_mat("HCPD", bquote(paste(italic("r"), " = ", .(round(HCPD_tractlevel_SA_ttest_p$rho.emp, 2)), ", ",
italic("p"[spin]), " = ", .(round(HCPD_tractlevel_SA_ttest_p$p.perm.cor, 3)))), binary = TRUE)
SA_age_mat_binary_HBN <- plot_meanSA_by_age_mat("HBN", bquote(paste(italic("r"), " = ", .(round(HBN_tractlevel_SA_ttest_p$rho.emp, 2)), ", ",
italic("p"[spin]), " = ", "N.S.")), binary = TRUE)
SA_age_mat_binary_avgdatasets <- plot_meanSA_by_age_mat("avgdatasets", bquote(paste(italic("r"), " = ", .(round(avgdatasets_tractlevel_SA_ttest_p$rho.emp, 2)), ", ", italic("p"[spin]), " < 0.0001")), binary = TRUE) + ggtitle("Average Across Datasets")
SA_age_mat_binary_plot_final <- ggarrange(SA_age_mat_binary_PNC, SA_age_mat_binary_HCPD, SA_age_mat_binary_HBN, SA_age_mat_binary_avgdatasets, ncol = 4, common.legend = T, legend = "bottom")
x.grob <- textGrob("Mean S-A Axis Rank of Cortical Endpoint",
gp=gpar(col="black", fontsize=24))
y.grob <- textGrob("Age of Maturation (Years)",
gp=gpar(col="black", fontsize=24), rot=90)
grid.arrange(arrangeGrob(SA_age_mat_binary_plot_final, left = y.grob, bottom = x.grob))aggregated_axis_PNC <- merge_SA_parcel("PNC") # 360 glasser parcels with mean age maturation and SA rank
aggregated_axis_HCPD <- merge_SA_parcel("HCPD")
aggregated_axis_HBN <- merge_SA_parcel("HBN")
# binarize age of maturation df's to get matured and not-yet-matured regions
# make an average binarized df (average across datasets first then binarize)
lh_PNC_merge <- lh_by_region_PNC %>% rename(PNC_mean_ageeffect = regional_mean_ageeffect)
lh_HCPD_merge <- lh_by_region_HCPD %>% rename(HCPD_mean_ageeffect = regional_mean_ageeffect)
lh_HBN_merge <- lh_by_region_HBN %>% rename(HBN_mean_ageeffect = regional_mean_ageeffect)
merge_temp <- merge(lh_HCPD_merge, lh_HBN_merge, by = "region")
lh_all_datasets_ageMat <- merge(merge_temp, lh_PNC_merge, by = "region")
lh_all_datasets_ageMat <- lh_all_datasets_ageMat %>% # left hemi age of maturation maps
mutate(regional_mean_ageeffect = rowMeans(select(., HCPD_mean_ageeffect, HBN_mean_ageeffect, PNC_mean_ageeffect), na.rm = TRUE))
rh_PNC_merge <- rh_by_region_PNC %>% rename(PNC_mean_ageeffect = regional_mean_ageeffect)
rh_HCPD_merge <- rh_by_region_HCPD %>% rename(HCPD_mean_ageeffect = regional_mean_ageeffect)
rh_HBN_merge <- rh_by_region_HBN %>% rename(HBN_mean_ageeffect = regional_mean_ageeffect)
merge_temp <- merge(rh_HCPD_merge, rh_HBN_merge, by = "region")
rh_all_datasets_ageMat <- merge(merge_temp, rh_PNC_merge, by = "region")
rh_all_datasets_ageMat <- rh_all_datasets_ageMat %>% # right hemi age of maturation maps
mutate(regional_mean_ageeffect = rowMeans(select(., HCPD_mean_ageeffect, HBN_mean_ageeffect, PNC_mean_ageeffect), na.rm = TRUE))
lh_all_datasets_ageMat$region <- paste0(lh_all_datasets_ageMat$region, "_L")
rh_all_datasets_ageMat$region <- paste0(rh_all_datasets_ageMat$region, "_R")
combined_df <- rbind(lh_all_datasets_ageMat, rh_all_datasets_ageMat) %>% rename(regionName = region)
binarized_combined_df <- right_join(combined_df, glasser_SAaxis, by = "regionName") # binarize
max_value <- max(binarized_combined_df$regional_mean_ageeffect, na.rm = TRUE)
binarized_combined_df$regional_mean_ageeffect[binarized_combined_df$regional_mean_ageeffect == max_value] <- NA
aggregated_axis_avgdatasets <- binarized_combined_df
# load pvalues
PNC_parcellevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/suppfigs_spintests/supp_parcellevel_pearson_pspin.csv")
HCPD_parcellevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/suppfigs_spintests/supp_parcellevel_pearson_pspin.csv")
HBN_parcellevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/suppfigs_spintests/supp_parcellevel_pearson_pspin.csv")
avgdatasets_parcellevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/avgdatasets/tract_profiles/suppfigs_spintests/supp_parcellevel_pearson_pspin_avgdatasets.csv")
# remember, this data is using bin = 5 (5 nodes at each end, after clipping 5 nodes already)
annot_text_PNC <- bquote(paste(italic("r"), " = ", .(round(PNC_parcellevel_SA_ttest_p$rho.emp, 2)), ", ",
italic("p"[spin]), " < 0.0001"))
annot_text_HCPD <- bquote(paste(italic("r"), " = ", .(round(HCPD_parcellevel_SA_ttest_p$rho.emp, 2)), ", ",
italic("p"[spin]), " < 0.0001"))
annot_text_HBN <- bquote(paste(italic("r"), " = ", .(round(HBN_parcellevel_SA_ttest_p$rho.emp, 2)), ", ",
italic("p"[spin]), " = ", .(round(HBN_parcellevel_SA_ttest_p$p.perm, 2))))
annot_text_avgdatasets <- bquote(paste(italic("r"), " = ", .(round(avgdatasets_parcellevel_SA_ttest_p$rho.emp, 2)), ", ",
italic("p"[spin]), " < 0.0001"))
agg_SA_parcel_PNC <- plot_agg_SA_parcel("PNC", annot_text_PNC, all_parcels = TRUE)
agg_SA_parcel_HCPD <- plot_agg_SA_parcel("HCPD", annot_text_HCPD, all_parcels = TRUE)
agg_SA_parcel_HBN <- plot_agg_SA_parcel("HBN", annot_text_HBN, all_parcels = TRUE)
agg_SA_parcel_avgdatasets <- plot_agg_SA_parcel("avgdatasets", annot_text_avgdatasets, all_parcels = TRUE) + ggtitle("Average Across Datasets")
agg_SA_parcel_plot_final <- ggarrange(agg_SA_parcel_PNC, agg_SA_parcel_HCPD, agg_SA_parcel_HBN, agg_SA_parcel_avgdatasets, ncol = 4)
x.grob <- textGrob("S-A Axis Rank", gp=gpar(col="black", fontsize=20))
y.grob <- textGrob("Age of Maturation (Years)", gp=gpar(col="black", fontsize=20), rot=90)
grid.arrange(arrangeGrob(agg_SA_parcel_plot_final, left = y.grob, bottom = x.grob))scalar = "dti_md"
PNC_ageeffects <- read.csv(paste0(output_root, "/", "PNC", "/GAM/", scalar, "/PNC_GAM_dev_measures.csv"))
HCPD_ageeffects <- read.csv(paste0(output_root, "/", "HCPD", "/GAM/", scalar, "/HCPD_GAM_dev_measures.csv"))
HBN_ageeffects <- read.csv(paste0(output_root, "/", "HBN", "/GAM/", scalar, "/HBN_GAM_dev_measures_withACT.csv"))# make a list of the different df names: [dataset]_GAM_ageeffects_[scalar]
datasets <- c("PNC", "HCPD", "HBN")
ageeffect_df_names <- c()
for (dataset in datasets) {
ageeffect_df_names <- append(ageeffect_df_names, paste0(dataset, "_ageeffects"))
}
# format age effect df's
ageeffect_dfs <- lapply(ageeffect_df_names, format_ageeffect)## There are 2273/2400 significant nodes
## There are 2392/2400 significant nodes
## There are 2175/2200 significant nodes
names(ageeffect.fdr_dfs) <- ageeffect_df_names
tracts <- setdiff(unique(ageeffect.fdr_dfs$HCPD_ageeffects$tract_label), "Corticospinal") # PNC
arc_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Arcuate_bin5_clip5_1sided.RData")
ifo_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Inferior_Fronto_occipital_bin5_clip5_1sided.RData")
ilf_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Inferior_Longitudinal_bin5_clip5_1sided.RData")
parc_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Posterior_Arcuate_bin5_clip5_1sided.RData")
slf_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Superior_Longitudinal_bin5_clip5_1sided.RData")
vof_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Vertical_Occipital_bin5_clip5_1sided.RData")
cc_af_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Anterior_Frontal_bin5_clip5_1sided.RData")
cc_mot_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Motor_bin5_clip5_1sided.RData")
cc_occ_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Occipital_bin5_clip5_1sided.RData")
cc_orb_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Orbital_bin5_clip5_1sided.RData")
cc_pp_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Posterior_Parietal_bin5_clip5_1sided.RData")
cc_sf_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Superior_Frontal_bin5_clip5_1sided.RData")
cc_sp_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Superior_Parietal_bin5_clip5_1sided.RData")
cc_tm_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Temporal_bin5_clip5_1sided.RData")
NEST_PNC <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC$pval, cc_mot_PNC$pval, cc_occ_PNC$pval, cc_orb_PNC$pval, cc_pp_PNC$pval, cc_sf_PNC$pval, cc_sp_PNC$pval, cc_tm_PNC$pval, arc_PNC$pval, ifo_PNC$pval, ilf_PNC$pval, parc_PNC$pval, slf_PNC$pval, vof_PNC$pval)), method=c("fdr")), rep("PNC", 14))) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))
# HCPD
arc_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Arcuate_bin5_clip5_1sided.RData")
ifo_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Inferior_Fronto_occipital_bin5_clip5_1sided.RData")
ilf_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Inferior_Longitudinal_bin5_clip5_1sided.RData")
parc_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Posterior_Arcuate_bin5_clip5_1sided.RData")
slf_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Superior_Longitudinal_bin5_clip5_1sided.RData")
vof_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Vertical_Occipital_bin5_clip5_1sided.RData")
cc_af_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Anterior_Frontal_bin5_clip5_1sided.RData")
cc_mot_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Motor_bin5_clip5_1sided.RData")
cc_occ_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Occipital_bin5_clip5_1sided.RData")
cc_orb_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Orbital_bin5_clip5_1sided.RData")
cc_pp_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Posterior_Parietal_bin5_clip5_1sided.RData")
cc_sf_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Superior_Frontal_bin5_clip5_1sided.RData")
cc_sp_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Superior_Parietal_bin5_clip5_1sided.RData")
cc_tm_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Temporal_bin5_clip5_1sided.RData")
NEST_HCPD <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_HCPD$pval, cc_mot_HCPD$pval, cc_occ_HCPD$pval, cc_orb_HCPD$pval, cc_pp_HCPD$pval, cc_sf_HCPD$pval, cc_sp_HCPD$pval, cc_tm_HCPD$pval, arc_HCPD$pval, ifo_HCPD$pval, ilf_HCPD$pval, parc_HCPD$pval, slf_HCPD$pval, vof_HCPD$pval)), method=c("fdr")), rep("HCP-D", 14))) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))
# HBN
arc_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Arcuate_bin5_clip5_1sided.RData")
ifo_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Inferior_Fronto_occipital_bin5_clip5_1sided.RData")
ilf_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Inferior_Longitudinal_bin5_clip5_1sided.RData")
parc_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Posterior_Arcuate_bin5_clip5_1sided.RData")
slf_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Superior_Longitudinal_bin5_clip5_1sided.RData")
vof_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Vertical_Occipital_bin5_clip5_1sided.RData")
cc_af_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Anterior_Frontal_bin5_clip5_1sided.RData")
cc_mot_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Motor_bin5_clip5_1sided.RData")
cc_occ_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Occipital_bin5_clip5_1sided.RData")
cc_orb_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Orbital_bin5_clip5_1sided.RData")
cc_pp_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Posterior_Parietal_bin5_clip5_1sided.RData")
cc_sf_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Superior_Frontal_bin5_clip5_1sided.RData")
cc_sp_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Superior_Parietal_bin5_clip5_1sided.RData")
cc_tm_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Temporal_bin5_clip5_1sided.RData")
NEST_HBN <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_HBN$pval, cc_mot_HBN$pval, cc_occ_HBN$pval, cc_orb_HBN$pval, cc_pp_HBN$pval, cc_sf_HBN$pval, cc_sp_HBN$pval, cc_tm_HBN$pval, arc_HBN$pval, ifo_HBN$pval, ilf_HBN$pval, parc_HBN$pval, slf_HBN$pval, vof_HBN$pval)), method=c("fdr")), rep("HBN", 14))) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))ageeffect_measure = "GAM.smooth.AdjRsq"
ageeffects_plots <- lapply(tracts, plot_age_effect, scalar = "dti_md", ageeffect_measure = ageeffect_measure,
color_HCPD = color_HCPD, color_HBN = color_HBN, color_PNC = color_PNC, clipEnds = 5, ylim2 = 0.41, ylim1 = 0, fontsize = 20, legend_position = "none")
names(ageeffects_plots) <- tractscolors <- c("#3FB8AFFF", "#FF9E9DFF", "#cc0468")
names(colors) <- c("HCP-D", "HBN", "PNC")
bin_num_nodes = 5
clipEnds = 5
# prepare the data for making lollipop plots
df_all <- bind_rows(
ageeffect.fdr_dfs$HCPD_ageeffects %>% mutate(Dataset = "HCP-D"),
ageeffect.fdr_dfs$HBN_ageeffects %>% mutate(Dataset = "HBN"),
ageeffect.fdr_dfs$PNC_ageeffects %>% mutate(Dataset = "PNC")
)
df_all$Dataset <- factor(df_all$Dataset, levels = c("HCP-D", "PNC", "HBN"))
df_formatted <- format_deep_superficial(df_all, ageeffect_measure, bin_num_nodes, clipEnds)
df_formatted <- df_formatted[complete.cases(df_formatted), ]
# merge with NEST results
NEST <- rbind(NEST_HCPD, NEST_HBN, NEST_PNC)
lollipop_data <- prepare_lollipop_data(df_formatted)
lollipop_data$Dataset <- factor(lollipop_data$Dataset, levels = c("PNC", "HCP-D", "HBN"))
# plot lollipops
tracts <- setdiff(unique(df_formatted$tract_label), "Corticospinal") # CST for supplement
lollipop_plots <- lapply(tracts, make_lollipop_plot, lollipop_data)
names(lollipop_plots) <- tracts# load my colormap :)
colormap_file <- sprintf("%1$s/code/colormaps/aquamarine.json", proj_root)
colormap <- fromJSON(paste(readLines(colormap_file), collapse=""))
aquamarine <- colorRampPalette(colormap)(15)
# arcuate
lh_arc_0.1 <- plot_threshold_maps(lh_maps_PNC$LeftArcuate, 'left', 0.1)
rh_arc_0.1 <- plot_threshold_maps(rh_maps_PNC$RightArcuate, 'right', 0.1)
arc_0.1 <- ggarrange(lh_arc_0.1, rh_arc_0.1, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
arc_0.1 <- annotate_figure(arc_0.1, top = text_grob("Arcuate Fasciculus", color = "black", size = 20))
lh_arc_0.3 <- plot_threshold_maps(lh_maps_PNC$LeftArcuate, 'left', 0.3)
rh_arc_0.3 <- plot_threshold_maps(rh_maps_PNC$RightArcuate, 'right', 0.3)
arc_0.3 <- ggarrange(lh_arc_0.3, rh_arc_0.3, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
arc_0.3 <- annotate_figure(arc_0.3, top = text_grob("Arcuate Fasciculus", color = "black", size = 20))
lh_arc_0.5 <- plot_threshold_maps(lh_maps_PNC$LeftArcuate, 'left', 0.5)
rh_arc_0.5 <- plot_threshold_maps(rh_maps_PNC$RightArcuate, 'right', 0.5)
arc_0.5 <- ggarrange(lh_arc_0.5, rh_arc_0.5, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
arc_0.5 <- annotate_figure(arc_0.5, top = text_grob("Arcuate Fasciculus", color = "black", size = 20))
# cst
lh_cst_0.1 <- plot_threshold_maps(lh_maps_PNC$LeftCorticospinal, 'left', 0.1)
rh_cst_0.1 <- plot_threshold_maps(rh_maps_PNC$RightCorticospinal, 'right', 0.1)
cst_0.1 <- ggarrange(lh_cst_0.1, rh_cst_0.1, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
cst_0.1 <- annotate_figure(cst_0.1, top = text_grob("Corticospinal Tract", color = "black", size = 20))
lh_cst_0.3 <- plot_threshold_maps(lh_maps_PNC$LeftCorticospinal, 'left', 0.3)
rh_cst_0.3 <- plot_threshold_maps(rh_maps_PNC$RightCorticospinal, 'right', 0.3)
cst_0.3 <- ggarrange(lh_cst_0.3, rh_cst_0.3, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
cst_0.3 <- annotate_figure(cst_0.3, top = text_grob("Corticospinal Tract", color = "black", size = 20))
lh_cst_0.5 <- plot_threshold_maps(lh_maps_PNC$LeftCorticospinal, 'left', 0.5)
rh_cst_0.5 <- plot_threshold_maps(rh_maps_PNC$RightCorticospinal, 'right', 0.5)
cst_0.5 <- ggarrange(lh_cst_0.5, rh_cst_0.5, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
cst_0.5 <- annotate_figure(cst_0.5, top = text_grob("Corticospinal Tract", color = "black", size = 20))
threshold_0.1 <- annotate_figure(ggarrange(arc_0.1 + theme(plot.margin = unit(c(0, 0, 1, 0), "cm")),
cst_0.1 + theme(plot.margin = unit(c(1, 0, 0, 0), "cm")),
nrow = 2, ncol = 1), top = text_grob("Threshold = 10%", color = "black", size = 20, face = 'bold'))
threshold_0.3 <- annotate_figure(ggarrange(arc_0.3 + theme(plot.margin = unit(c(0, 0, 1, 0), "cm")),
cst_0.3 + theme(plot.margin = unit(c(1, 0, 0, 0), "cm")),
nrow = 2, ncol = 1), top = text_grob("Threshold = 30%", color = "black", size = 20, face = 'bold'))
threshold_0.5 <- annotate_figure(ggarrange(arc_0.5 + theme(plot.margin = unit(c(0, 0, 1, 0), "cm")),
cst_0.5 + theme(plot.margin = unit(c(1, 0, 0, 0), "cm")),
nrow = 2, ncol = 1), top = text_grob("Threshold = 50%", color = "black", size = 20, face = 'bold'))
all_threshold_plots <- ggarrange(threshold_0.1, threshold_0.3, threshold_0.5, ncol = 3, nrow = 1)
all_threshold_plotsPNC_lh_GM <- read.csv("/cbica/projects/luo_wm_dev/input/PNC/derivatives/gyral_hops/all_subjects/lh_GM_count.csv")
PNC_rh_GM <- read.csv("/cbica/projects/luo_wm_dev/input/PNC/derivatives/gyral_hops/all_subjects/rh_GM_count.csv")
PNC_lh_WM <- read.csv("/cbica/projects/luo_wm_dev/input/PNC/derivatives/gyral_hops/all_subjects/lh_WM_count.csv")
PNC_rh_WM <- read.csv("/cbica/projects/luo_wm_dev/input/PNC/derivatives/gyral_hops/all_subjects/rh_WM_count.csv")
# 1.5 mm is best
kbl(PNC_lh_GM, align = "lccrr", caption = "PNC Left Hemi: best depths for GM") %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "left") | Depth | count |
|---|---|
| depth_0 | 0 |
| depth_0p5 | 0 |
| depth_1 | 375 |
| depth_1p5 | 726 |
| depth_2 | 0 |
kbl(PNC_rh_GM, align = "lccrr", caption = "PNC Right Hemi: best depths for GM") %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "left") | Depth | count |
|---|---|
| depth_0 | 0 |
| depth_0p5 | 0 |
| depth_1 | 385 |
| depth_1p5 | 716 |
| depth_2 | 0 |
kbl(PNC_lh_WM, align = "lccrr", caption = "PNC Left Hemi: best depths for WM") %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "left") | Depth | count |
|---|---|
| depth_0 | 0 |
| depth_0p5 | 0 |
| depth_1 | 272 |
| depth_1p5 | 829 |
| depth_2 | 0 |
kbl(PNC_rh_WM, align = "lccrr", caption = "PNC Right Hemi: best depths for WM") %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "left") | Depth | count |
|---|---|
| depth_0 | 0 |
| depth_0p5 | 0 |
| depth_1 | 273 |
| depth_1p5 | 828 |
| depth_2 | 0 |
tracts <- c("Callosum Anterior Frontal", "Callosum Motor", "Callosum Occipital", "Callosump Orbital", "Callosum Posterior Parietal", "Callosum Superior Frontal", "Callosum Superior Parietal", "Callosum Temporal", "Arcuate", "Inferior Fronto-occipital", "Inferior Longitudinal", "Posterior Arcuate", "Superior Longitudinal", "Uncinate", "Vertical Occipital")
# PNC
# bin size = 3
arc_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Arcuate_bin3_clip5_1sided.RData")
ifo_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Fronto_occipital_bin3_clip5_1sided.RData")
ilf_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Longitudinal_bin3_clip5_1sided.RData")
parc_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Posterior_Arcuate_bin3_clip5_1sided.RData")
slf_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Superior_Longitudinal_bin3_clip5_1sided.RData")
unc_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Uncinate_bin3_clip5_1sided.RData")
vof_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Vertical_Occipital_bin3_clip5_1sided.RData")
cc_af_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Anterior_Frontal_bin3_clip5_1sided.RData")
cc_mot_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Motor_bin3_clip5_1sided.RData")
cc_occ_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Occipital_bin3_clip5_1sided.RData")
cc_orb_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Orbital_bin3_clip5_1sided.RData")
cc_pp_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Posterior_Parietal_bin3_clip5_1sided.RData")
cc_sf_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Frontal_bin3_clip5_1sided.RData")
cc_sp_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Parietal_bin3_clip5_1sided.RData")
cc_tm_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Temporal_bin3_clip5_1sided.RData")
NEST_PNC_bin3 <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC_bin3$pval, cc_mot_PNC_bin3$pval, cc_occ_PNC_bin3$pval, cc_orb_PNC_bin3$pval, cc_pp_PNC_bin3$pval, cc_sf_PNC_bin3$pval, cc_sp_PNC_bin3$pval, cc_tm_PNC_bin3$pval, arc_PNC_bin3$pval, ifo_PNC_bin3$pval, ilf_PNC_bin3$pval, parc_PNC_bin3$pval, slf_PNC_bin3$pval, unc_PNC_bin3$pval, vof_PNC_bin3$pval)), method=c("fdr")), rep("PNC", 15))) %>% setNames(c("tract_label", "bin_3_pvalue", "Dataset"))
# bin size = 5
arc_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Arcuate_bin5_clip5_1sided.RData")
ifo_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Fronto_occipital_bin5_clip5_1sided.RData")
ilf_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Longitudinal_bin5_clip5_1sided.RData")
parc_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Posterior_Arcuate_bin5_clip5_1sided.RData")
slf_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Superior_Longitudinal_bin5_clip5_1sided.RData")
unc_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Uncinate_bin5_clip5_1sided.RData")
vof_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Vertical_Occipital_bin5_clip5_1sided.RData")
cc_af_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Anterior_Frontal_bin5_clip5_1sided.RData")
cc_mot_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Motor_bin5_clip5_1sided.RData")
cc_occ_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Occipital_bin5_clip5_1sided.RData")
cc_orb_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Orbital_bin5_clip5_1sided.RData")
cc_pp_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Posterior_Parietal_bin5_clip5_1sided.RData")
cc_sf_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Frontal_bin5_clip5_1sided.RData")
cc_sp_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Parietal_bin5_clip5_1sided.RData")
cc_tm_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Temporal_bin5_clip5_1sided.RData")
NEST_PNC_bin5 <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC_bin5$pval, cc_mot_PNC_bin5$pval, cc_occ_PNC_bin5$pval, cc_orb_PNC_bin5$pval, cc_pp_PNC_bin5$pval, cc_sf_PNC_bin5$pval, cc_sp_PNC_bin5$pval, cc_tm_PNC_bin5$pval, arc_PNC_bin5$pval, ifo_PNC_bin5$pval, ilf_PNC_bin5$pval, parc_PNC_bin5$pval, slf_PNC_bin5$pval, unc_PNC_bin5$pval, vof_PNC_bin5$pval)), method=c("fdr")), rep("PNC", 15))) %>% setNames(c("tract_label", "bin_5_pvalue", "Dataset"))
# bin size = 7
arc_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Arcuate_bin7_clip5_1sided.RData")
ifo_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Fronto_occipital_bin7_clip5_1sided.RData")
ilf_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Longitudinal_bin7_clip5_1sided.RData")
parc_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Posterior_Arcuate_bin7_clip5_1sided.RData")
slf_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Superior_Longitudinal_bin7_clip5_1sided.RData")
unc_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Uncinate_bin7_clip5_1sided.RData")
vof_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Vertical_Occipital_bin7_clip5_1sided.RData")
cc_af_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Anterior_Frontal_bin7_clip5_1sided.RData")
cc_mot_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Motor_bin7_clip5_1sided.RData")
cc_occ_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Occipital_bin7_clip5_1sided.RData")
cc_orb_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Orbital_bin7_clip5_1sided.RData")
cc_pp_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Posterior_Parietal_bin7_clip5_1sided.RData")
cc_sf_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Frontal_bin7_clip5_1sided.RData")
cc_sp_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Parietal_bin7_clip5_1sided.RData")
cc_tm_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Temporal_bin7_clip5_1sided.RData")
NEST_PNC_bin7 <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC_bin7$pval, cc_mot_PNC_bin7$pval, cc_occ_PNC_bin7$pval, cc_orb_PNC_bin7$pval, cc_pp_PNC_bin7$pval, cc_sf_PNC_bin7$pval, cc_sp_PNC_bin7$pval, cc_tm_PNC_bin7$pval, arc_PNC_bin7$pval, ifo_PNC_bin7$pval, ilf_PNC_bin7$pval, parc_PNC_bin7$pval, slf_PNC_bin7$pval, unc_PNC_bin7$pval, vof_PNC_bin7$pval)), method=c("fdr")), rep("PNC", 15))) %>% setNames(c("tract_label", "bin_7_pvalue", "Dataset"))
# bin size = 10
arc_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Arcuate_bin10_clip5_1sided.RData")
ifo_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Fronto_occipital_bin10_clip5_1sided.RData")
ilf_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Longitudinal_bin10_clip5_1sided.RData")
parc_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Posterior_Arcuate_bin10_clip5_1sided.RData")
slf_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Superior_Longitudinal_bin10_clip5_1sided.RData")
unc_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Uncinate_bin10_clip5_1sided.RData")
vof_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Vertical_Occipital_bin10_clip5_1sided.RData")
cc_af_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Anterior_Frontal_bin10_clip5_1sided.RData")
cc_mot_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Motor_bin10_clip5_1sided.RData")
cc_occ_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Occipital_bin10_clip5_1sided.RData")
cc_orb_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Orbital_bin10_clip5_1sided.RData")
cc_pp_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Posterior_Parietal_bin10_clip5_1sided.RData")
cc_sf_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Frontal_bin10_clip5_1sided.RData")
cc_sp_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Parietal_bin10_clip5_1sided.RData")
cc_tm_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Temporal_bin10_clip5_1sided.RData")
NEST_PNC_bin10 <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC_bin10$pval, cc_mot_PNC_bin10$pval, cc_occ_PNC_bin10$pval, cc_orb_PNC_bin10$pval, cc_pp_PNC_bin10$pval, cc_sf_PNC_bin10$pval, cc_sp_PNC_bin10$pval, cc_tm_PNC_bin10$pval, arc_PNC_bin10$pval, ifo_PNC_bin10$pval, ilf_PNC_bin10$pval, parc_PNC_bin10$pval, slf_PNC_bin10$pval, unc_PNC_bin10$pval, vof_PNC_bin10$pval)), method=c("fdr")), rep("PNC", 15))) %>% setNames(c("tract_label", "bin_10_pvalue", "Dataset"))PNC_bin_compare <- left_join(NEST_PNC_bin3, NEST_PNC_bin5) %>% left_join(NEST_PNC_bin7) %>% left_join(NEST_PNC_bin10) %>% select(-Dataset)
PNC_bin_compare$bin_3_pvalue <- as.numeric(PNC_bin_compare$bin_3_pvalue)
PNC_bin_compare$bin_5_pvalue <- as.numeric(PNC_bin_compare$bin_5_pvalue)
PNC_bin_compare$bin_7_pvalue <- as.numeric(PNC_bin_compare$bin_7_pvalue)
PNC_bin_compare$bin_10_pvalue <- as.numeric(PNC_bin_compare$bin_10_pvalue)
# binning differences don’t impact significance status, though specific p-value changes a little
PNC_bin_table <- PNC_bin_compare %>%
mutate(across(where(is.numeric), ~ formatC(.x, digits = 2, format = "fg", drop0trailing = TRUE))) %>%
rename(Tract = tract_label,
"Bin Size 3" = bin_3_pvalue,
"Bin Size 5" = bin_5_pvalue,
"Bin Size 7" = bin_7_pvalue,
"Bin Size 10" = bin_10_pvalue)
PNC_bin_table_final <- PNC_bin_table %>%
kbl(format = "html", align = "lcccc", digits = 2) %>%
add_header_above(c(" " = 1, "p-value" = 4), bold = TRUE) %>%
kable_styling(full_width = FALSE, position = "center") %>%
kable_classic(full_width = FALSE, html_font = "Arial") %>%
row_spec(0:nrow(PNC_bin_table), extra_css = "line-height: 0.5;")
PNC_bin_table_final| Tract | Bin Size 3 | Bin Size 5 | Bin Size 7 | Bin Size 10 |
|---|---|---|---|---|
| Callosum Anterior Frontal | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Callosum Motor | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Callosum Occipital | 0.00012 | 0.00012 | 0.00012 | 0.032 |
| Callosump Orbital | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Callosum Posterior Parietal | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Callosum Superior Frontal | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Callosum Superior Parietal | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Callosum Temporal | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Arcuate | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Inferior Fronto-occipital | 0.00012 | 0.00012 | 0.0093 | 0.0075 |
| Inferior Longitudinal | 0.36 | 0.22 | 0.26 | 0.24 |
| Posterior Arcuate | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Superior Longitudinal | 0.00012 | 0.00012 | 0.00012 | 0.00014 |
| Uncinate | 0.17 | 0.15 | 0.12 | 0.1 |
| Vertical Occipital | 0.00012 | 0.00012 | 0.00012 | 0.00014 |